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ABSTRACT 


An  efficient  .method  is  given  for  computing  the  incomplete  beta 
function  ratio,  Ix(a,b),  on  a  high  speed  digital  computer.  The 
arguments  a,b,  are  limited  to  positive  integral  multiples  of  one- 
half  values  over  the  ranges  l/2  ^  a  ^  10s  ,  l/2  ^  b  ^  60. 

The  program  has  been  coded  in  STRAP  for  the  IBM  7030  (STRETCH). 
The  average  computing  time  for  a  ten  decimal  digit  value  of  Ix(a,b) 
is  2.6  milliseconds;  on  an  IBM  7090  it  would  be  about  8  milliseconds 
per  case. 
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FOREWORD 


The  work  reported  in  this  publication  was  done  in  the  Applied 
Mathematics  Section  of  the  Mathematics  Research  Group  with  Foundational 
Research  funds  No.  R360FR103/2101/R0110101 . 

The  development  of  an  Ix(a,b)  routine  was  requested  by  Dr.  K.  Abt 
of  the  Mathematical  Statistics  Branch,  Operations  Research  Division. 

The  routine  is  of  vital  importance  as  a  subroutine  in  a  larger  sta¬ 
tistical  program  (NOVACOM) ,  presently  under  development.  NOVAGOM, 
which  will  perform  analysis  of  variance  for  data  classifications  with 
missing  observations,  is  a  program  of  wide  applicability  in  weapons 
effectiveness  studies  and  other  statistical  problems. 

The  IBM  7030  code  was  developed  by  Mr.  Travis  Herring  from  flow 
charts  contained  herein.  Auxiliary  subroutines  for  the  calculation  of 
certain  elementary  functions  were  taken  from  the  7030  Systems  Library 
subroutines . 

A  NORC  code  for  ^(ajb)  was  initially  constructed  primarily  for 
exploratory  type  calculations.  The  auxiliary  subroutines  incorporated 
in  the  NORC  code  were  taken  from  the  library  of  NORC  subroutines 
developed  by  Dr.  A.  V.  Hershey  of  the  Science  Research  Group. 

This  report  contains  more  recent  developments  and  supersedes 
NWL  Report  1949  of  28  February  1965. 

Date  of  completion  was  October  1966. 


APPROVED  FOR  RELEASE: 


/s/  RALPH  A.  NIEMANN,  Acting 
Technical  Director 
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I.  INTRODUCTION 

The  incomplete  beta  function,  Bx(a,b)  is  defined  as  follows: 


x 

Bx(a,b)  E  J  ta 

O 

where 

0  ^  x  £  1  , 

For  x  =  1,  Bx(a,b)  is  known  as 
be  expressed  in  terms  of  three 


1  (1  -  t)10-1  dt  ,  (1) 


a  >  0  ,  b>0. 

the  complete  beta  function.  It  can 
complete  gamma  functions,  [2;  p  127], 


B  (a,b)  =  ,  (2) 

1  r(a  +  b) 

where  the  complete  gamma  function,  with  argument  s. ,  is  defined  by 


P(s) 


s  >  0  . 


(3) 


Throughout  this  paper  the  following  constraints  are  imposed  on  a 
and  b: 


{1}  They  can  only  tahe  positive  values  of  integral  multiples 
of  one -half, 

{2}  They  satisfy  the  inequalities: 

l/2  <  a  s  108  ,  l/2  <  b  <;  60  . 

The  second  constraint  may  be  relaxed  on  the  upper  bound  of  b,  at  a 
proportional  increase  in  the  amount  of  calculation  required. 

The  purpose  of  this  report  is  to  describe  an  efficient 
method  on  a  digital  computer  for  the  high  accuracy  computation  of 
the  ratio  of  (l)  to  (2)  subject  to  the  constraints  {1}  and  {2)  on 


1 


ML  REPORT  NO.  1949 


a  and  1.  This  ratio  is  known  as  the  incomplete  beta  function 
ratio  and  is  indicated  by  the  symbolism  Ix(a,b),  i.e., 

Ix(a,b)  =  Bx(a,b)/B1(a,b)  .  (4) 

By  the  substitution  of  u  =  1  -  t  in  (l) 

Ix(a,b)  =  1  -  I1-X(b;a)  .  (5) 

In  probability  terms,  Ix(a,b)  is  called  the  beta  distribu¬ 
tion  function,  with  mean  p  and  variance  a2  given  by 

p  =  a/(a  +  b)  ,  a2  =  ab/[(a  +  b  +  l)(a  +  b)2]  ,  [2;  p  244].  (6) 

The  importance  of  this  function  is  reflected  in  Karl 
Pearson*  s  monumental  work.  Tables  of  the  Incomplete  Beta  Function, 
[10],  which  required  ten  years  to  complete  (1923-1932).  The 
method  he  employed  will  be  outlined  in  Section  III.  The  primary 
importance  of  the  beta  distribution  function,  Ix(a,b),  stems  from 
the  fact  that  it  is  directly  related  or  interpreted  in  terms  of 
three  basic  continuous  probability  distribution  functions,  the 
chi-square  distribution,  the  F  (variance  ratio)  distribution,  and 
the  Student *s  t  distribution.  It  is  also  related  to  the  discrete 
cumulative  binomial  distribution. 

It  will  be  shown  in  the  subsequent  discussion  that  the 
constraint  {1}  above  is  not  a  very  severe  one,  since  the  import¬ 
ant  related  distributions,  just  mentioned,  are  covered  by  the 
values  of  a  and  b  allowed  under  {1}. 

The  remainder  of  this  section  is,  for  the  most  part,  taken 
from  [i;  p  940-948].  Let  X±9  X2,  ...  Xv  be  independent  and 
identically  distributed  random  variables  each  following  a  normal 
distribution  with  mean  zero  and  unit  variance.  Then,  as  it  is 

v 

known  in  statistics,  X2  =  Z  X2  is  said  to  follow  a  chi-square 

i  =  i 

distribution  with  v  degrees  of  freedom;  the  probability  of  the 
event  X2  ^  X2  is  given  by 

X2 

P(X2  2  X2}  =  P(X2|  v)  =  [2 v^2  rCv/S)]'1  J  e”t/2  pvAO-i  dt  ,  (7) 

o 


A  proof  of  (7)  is  given  in  [2;  p  233]. 
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Now,  if  X2  and  x|  are  independent  random  variables  which 
follow  chi-square  distributions  with  vx  and  degrees  of  freedom 
respectively,  thbn  Xf/(X^  +  X2)  follows  a  beta  distribution  where 
a  =  Vj/2,  b  =  v2/2.  Thus 

P{X2/(X2  +xj)  sx}  =  Ix(a,b)  .  (8) 


A  proof  of  (8)  is  given  in  [2;  p  243]. 

If  we  consider  the  same  random  variables  X2,  X2,  then  the 
distribution  of  the  ratio 


F  = 


xf/vo 


(9) 


2^  *2 

is  said  to  follow  the  variance  ratio  or  F  distribution  with  v.]_ 
and  vP  degrees  of  freedom.  The  probability  that  F  s  F  is  given  by 


P{F  £  F0]  =P(F|v,,V2) 


o'  Yl>  v2> 
VW2 


V2  „  V2 


F0  /V 


Bx(  Vi/2,  V2/2) 


*  0  . 


-(v  +V.  )/2 

2  (v2  +  V:  F)  21  dF  , 


(10) 


A  proof  of  (10)  is  given  in  [2;  p  241-3],  The  substitution 
z  =  v2/(v2  +  v  F)  is  applied  to  (10).  It  follows  directly 

from  this  variable  of  integration  substitution  that 

f  Vi  Vo\ 


where 


x  = 


V2  +  vx  Fc 


(2k 

XV  ** 

t)  -  M 

\  2 

2  '  \ 

V 

> 

1 

11 

X 

H 

V2 

(11) 
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If  Xj_  is  a  random  variable  following  a  normal  distribution 

with  mean  zero  and  unit  variance,  and  X2  is  a  random  variable 
following  an  independent  chi-square  distribution  with  v  degrees 
of  freedom,  the  distribution  of  the  ratio 

xJ/vFJv  =  t  (12) 

is  called  the  Student's  t-distribution  with  v  degrees  of  freedom. 
The  probability  that  |t|  will  be  less  than  a  fixed  constant ,  t  , 
is  given  by 

p(  bl  *  t0)  =  A(tJ  v) 

r  /  n-i  p°  v+i 

=  '^Bi(|^|)  J  (1  +  t2/v)"  2  dt  .  (13) 


A  proof  of  (13)  is  given  in  [2;  p  237-240],  In  terms  of  the  beta 
distribution 


A(t0lv)  =  1  -  lii. 


V  +  t! 


The  derivation  is  straightforward;  apply  the  transformation 
z  =  v/(v  +  t2)  to  (13). 

In  case  a  and  b  are  specified  as  positive  integers,  the  beta 
function  is  related  directly  to  the  cumulative  binomial  distribu¬ 
tion,  E(n,r,x),  which  is  given  by 


n 

E(n,r,x)  =  Z  e(n,i,x) 
i=r 


(15) 


where 


e(n,i,x)  h  ^  x1(l  -  x)51-1;  [21],  [22.]  .  (16) 
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If  x  is  the  probability  of  success  in  one  trial,  the  cumulative 
binomial  distribution,  E(n,r,x),  represents  the  probability  that 
at  least  r  successes  will  occur  in  n  independent  trials,  and 
e(n,r,x)  the  probability  that  there  will  occur  exactly  r  successes 
in  n  independent  trials.  In  terms  of  Ix,  we  have 

B  (r,n  -  r  +  l)  .  .  ,  x 

E(n,r,x)  =  JL -  =  I  (r,n  -  r  +  l)  .  (17) 

B j( r , n  -  r  +  1) 

The  derivations  of  (15)  and  (16)  are  given  in  [2;  p  193-194];  (17) 
is  derived  in  [21;  p  XVII].  Applications  of  E(n,r,x)  are  also  given 
in  [21,  p  XXXIV];  applications  of  the  continuous  distributions 
given  above  can  be  found  in  many  references,  e.g.,  throughout  [5]. 

The  four  distributions  that  have  just  been  related  to  the 
beta  distribution  require  only  positive  values  of  a  and  b  at  the 
integral  multiples  of  one-half.  Moreover,  even  the  non- central  F 
and  t  distributions  are  included  by  these  values  of  a  and  b, 

[1;  p  947]. 

A  number  of  published  tables  exist  for  I  or  its  inverse  with 
respect  to  x.  Four  of  the  most  extensive  ones  are  referenced  here 
and  the  ranges  of  the  variables  are  given. 

A  table  already  mentioned  is  K.  Pearson' s, [10] .  It  is  the 
most  comprehensive  for  Ix(a,b).  The  ranges  of  a,b,  and  x  are 

a  =  ^^)50,  b  =  7^7^50  such  that  a  ^  b,  x  =  0(  .01)1.00.  Values  of 

the  beta  ratio  are  printed  to  seven  decimal  digits. 

For  integer  values  of  a  and  b  there  are  the  Tables  of  the 
Cumulative  Binomial  Probability  Distribution,  [21],  issued  by 
Harvard  University  in  1955,  and  the  Tables  of  the  Binomial  Prob¬ 
ability  Distribution,  [22  ],  published  by  the  Bureau  of  Standards  in 
1950.  The  ranges  of  n,r,x  in  [21]  are: 

r  =  0(l)n,  n  =  1(1)50(2)100(10)200(20)500(50)1000, 

x  =  0.01(0.01)0.50,  1/16,  1/12,  1/8,  l/6,  3/16,  5/16,  1/3,  3/8,  5/12, 
7/16. 
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Ix(r,n  -  r  +  1)  is  given  to  five  decimal  digits,  where  a  =  r, 
h  =  n  -  r  +  1.  The  ranges  of  n,r,x  in  [22]  are: 

x  =  0.01(0.01)0.50,  n  =  2(l)49,  r  =  l(l)n  , 

and  Ix(r,n  -  r  +  l)  is  given  to  seven  decimal  digits. 

A  table  of  percentage  points  of  Ix(a,b)  has  been  computed  by 
C.  M.  Thompson,  [15].  In  this  case,  the  variable  x  is  tabulated 
as  a  function  of  Ix,a,b.  The  ranges  are: 

Ix  =  0.005,  .01,  .025,  .05,  0.1,  .25,  .5;  2a  =  l(l)30,  40,  60,  120,  » 

2b  =  1(1)10,  12,  15,  20,  24,  30,  40,  60,  120  . 

The  computed  value  of  x  is  given  to  five  decimal  digits. 

In  [9]  a  nomogram  is  given  which  yields  graphical  results  for 
Ix(a,b)  somewhat  beyond  the  ranges  of  [10]  and  [15].  The  values  of 
a  and  b  extend  to  70  and  60  respectively. 
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The  importance  of  Ix(a,b)  makes  it  extremely  desirable  to  have 
a  digital  computer  program  which  is  designed  for  the  efficient  cal¬ 
culation  of  Ix(a,b)  to  say  eight  or  more  decimal  digits  for  any 

values  of  a  and  b  subject  to  the  constraints  on  page  1.  To  the 
authors'  knowledge,  no  efficient  program  exists  for  such  a  calcula¬ 
tion.  A  description  of  a  program  which  is  suitable  is  given  in 
Sections  IV  and  V. 

The  next  section  includes  a  discussion  of  some  previously  pub¬ 
lished  formulas,  algorithms  and  computing  programs.  In  order  to 
more  easily  set  forth  where  some  of  these  methods  fail  to  be  useful, 
the  major  numerical  difficulties  in  computing  Ix  are  stated: 

{a}  A  straightforward  binomial  expansion  of  the  integrand 
in  (l)  and  a  subsequent  integration  to  obtain  an  alternating  series 
in  powers  of  x  cannot  be  used  for  large  values  of  a  and  b.  The 
eventual  subtraction  of  consecutive  terms  of  nearly  equal  absolute 
value  causes  a  loss  in  significant  digits  -which  is  prohibitive. 

fb}  Ix(a,b)  is  a  function  of  three  independent  variables.  It 
is  unlikely  that  one  procedure  or  algorithm  will  suffice,  and  so  it 
will  be  necessary  to  devise  a  variety  of  techniques  over  the  ranges 
of  a,  b,  and  x  which  are  contemplated. 

(c)  The  extensive  range  of  a,  |  £  a  <:  10s,  introduces  scaling 

problems  in  most  procedures  because  terms  of  the  order  of  r(a) 
occur. 


fd)  The  use  of  recurrence  relations  imposes  the  requirement 
of  computing  starting  values,  in  which  case,  one  is  confronted  with 
the  evaluation  of  quantities  such  as  Ix(a,l/2)  for  large  a.  This 
computation  is  not  straightforward. 

{e}  Closely  connected  to  {d}  is  the  fact  that  one  must  dodge 
any  procedure  which  attempts  to  sum  over  a  elements,  since  this 
could  entail  the  addition  of  10s  elements.  Such  a  process  would 
destroy  the  efficiency  of  the  program  and  very  likely  the  accuracy 
as  well. 
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III.  VARIOUS  METHODS  FOR  COMPUTING  I  (a,b),  iji,  b 

x  xy2 


A  search  for  finding  a  suitable  set  of  methods  for  computing 
Ix(a,b)  was  initiated  by  carrying  out  an  investigation  of  the 

literature  on  the  subject.  In  this  section  some  of  the  more  perti¬ 
nent  papers,  from  our  point  of  view,  are  discussed,  and  reasons  for 
not  using  a  particular  program  or  analysis  are  pointed  out. 

It  seems  fitting  to  begin  with  the  algorithms  used  by  K. 
Pearson  in  computing  his  table,  [10].  They  are  founded  on  the 
recurrence  relations: 

Ix(a,b)  =  x  Ix(a  -  l,b)  +  (l  -  x)  Ix(a,b  -  l)  ,  (18) 


Ix(a  +  1,  1/2)  =  Ix(a,3/2) 


2r(a  +  5/2)  xa  /I  -  x 
v/tt  r(a  +  l) 


Ix(l/2,  b  +  1)  =  Ix(3/2,b)  +  gr(ft  +  3/2)  7x  (l  -  x) 

/rf  r(b  +  l) 


which  are  derived  in  Appendix  A. 

The  procedure  is  exempli¬ 
fied  graphically.  In  Figure  (l) 
the  order  of  the  computation  for 
the  case  of  integral  multiples 
of  one-half  for  a  and  b,  where 
a  and  b  extend  as  far  as  7/2,  is 
shown.  The  ordered  pair  of 
values  (a,b)  at  a  node  are  those 
for  which  Ix(a,b)  is  computed. 
The  circled  digits  specify  the 
order  in  which  the  consecutive 
values  of  Ix  are  obtained. 
Although  this  process  is  ade¬ 
quate  for  generating  a  table  of 
Ix,  it  certainly  would  not  be 
efficient  if  a  exceeded  50  by 
any  significant  amount.  Diffi¬ 
culties  fd)  and  [e]  of  the  last 
section  would  be  encountered. 
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Numerous  formulas  have  been  developed  and  investigated  by 
H.  E.  Soper  and  reported  in  [12].  His  work  includes  perhaps  a 
greater  variety  of  relations  for  I  than  any  other  published  paper 

on  the  subject.  The  contents  include  the  rates  of  convergence  of 
derived  series  for  arbitrary  real  positive  values  of  a  and  b. 
Polynomial  fits  and  Fourier  series  expansions  are  considered.  The 
general  conclusion  of  Soper  appears  to  be  however,  that  none  of 
the  methods  given,  other  than  those  used  by  K.  Pearson  are  ade¬ 
quate  for  computation,  [12;  p  49]. 


A  paper  by  J.  Wishart,  [20],  resolves  difficulty  [c]  for 
sufficiently  large  a  and  b.  For  completeness1  sake,  a  derivation 
of  his  results  is  given. 


By  (1) 

n 


Ml’* 


sin""1  ,/x 

(l  -  t)b_1  dt//t  =  2  f  (cos  e)2b_1  <3.0  ,  (21) 


where  t  =  sin2  6.  If  /t  =  2y/(l  +  y2) ,  then 


B, 


4  /  u . .  y h2b~V  ay 

2>  l  1  +  y2 


i  +  r 


X 


(22) 


r 


=  4  /  exp  j  -  2W  z  (y2)21-V(2i  -  l) 
J  L  1  =  1 


dy 


1  +  ys 


where 


=  exp 

■ 

-  R  In 

f1 '  Z)] 

Vl  +  y2/ 

\1  +  y2/. 
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If  y  =  z/(2  ,/TjJ  and  if  the  power  series  in  z2/4N  is  used  for 
£n[(l  -  y2)/(l  +  y2)],  then  (22)  becomes 


2X  /N 


Bxfi,  =  *.  r 

\2  )  JRJ  1 


1  +  (z/2N)2 
2X  /n 


00  /  2  \2i  -  1  ✓,  .11 

'  -  2N  Z  /  (21  ■  1)  dz 


\  4N  16N2 


z6  z10  .  z12  \  , 

- —  -  - -  +  - —  +  ...  dz 

96N2  2560 N4  18432R6  / 


2X  /n 


i  +  i) 

3  2N  / 


w«l 

1 

H 

t  ,, 

z4  . 

[  4N 

16N2 

z1Q 

z12 

960N4 

18432N4 

128R3  '3  2N  /  960N4  18432N4 J 

Equation  (24)  is  obtained  by  expanding  the  product  of  exponentials 
in  (23)  in  powers  of  z2.  The  approximation  (25)  is  the  result  of 
truncating  the  power  series  in  (24)  to  terms  of  order  l/N4.  The 

indicated  term  by  term  integration  of  (25)  and  division  by  b^ 

gives 


J4i’ 


717 L"c  -  45“=  <p) 


-^5  (-  +  — )  m6(p)  + 

32N2  V3  2N/ 


105  /I 


63  „  ,  1155 

7  mio(P)  +  " 

>4N4  2048N4 


128N3  \ 3  2N 
m1s>(  B)  ) 


M-S(P) 
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where  (3=2  /N  X  and  the  incomplete  normal  moment  function, 


P 

m  (p)  =  JL  f  z1  e~z  /2  dz/( i  -  l)  (i  -  3) 

'/2tTJo 


2  or  1  . 


(27) 


Regrouping  terms,  Wishart's  final  result  is  obtained. 


(28) 


where 

<J0(p)  =  m0(p),  d>1  =  (l/4)  ^(p);  <t>2( P)  *  *1875  m4(p)  “  *15625  mg(p)  , 
0  (p)  s  .234375  m6(p)  -  .2734375  mQ(p)  , 

04(p)  ■  .41015625  m8(p)  -  .984375  m10(p)  +  .56396484  m12(p)  , 


and 


8=2  /(2b  -  +  /l  -  x)  . 


(29) 


Equation  (28)  has  the  very  desirable  feature  of  approaching 


the  correct  limiting  value  for  I: 


as  b 


The  equation 


was  not  employed  though,  because  it  would  have  required  incorpora¬ 
tion  into  our  program  of  an  efficient  normal  probability  integral 
subroutine.  A  fast  subroutine  for  the  probability  integral  gener¬ 
ally  requires  storage  of  a  set  of  function  values  at  the  expense  of 
300  to  500  storage  locations  in  the  computer.  Also,  a  great  deal 
more  numerical  analysis  would  have  been  required  on  (28)  to  fix 
rigorous  error  bounds  and  to  determine  the  range  over  which  it  could 
be  used  efficiently.  It  is  difficult  to  decide  without  the  addi¬ 
tional  study  whether  it  would  be  worthwhile  to  insert  the  procedure 
into  our  present  program,  especially  since  the  procedure  we  employ 
for  this  calculation  is  quite  efficient  (see  Sections  IV.  and VI)  in 
its  own  right. 
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H.  E.  Fettis,  [3],  treats  the  problem  of  evaluating 


sinN  0  cosM  0  d0  numerically.  This  integral  corresponds  to 


o 

1  b  [N_±l  |  where  x  =  sin2  0.  It  is  emphasized  at  the 

2  x\  2  2  ) 

outset  that  it  is  believed  Fettis  was  only  interested  in  arbitrary 
small  positive  values  of  N  and  M.  Nevertheless  it  was  not  obvious 
at  first  sight  whether  his  formula  would  be  useful  for  large  N 
and  M.  As  it  turns  out,  they  are  not.  The  basic  formula  of  the 
paper  is  given  by 


0 


M 

0  cos 


0  d0 


0 


8  [l  .  (M.-  1  ?  (»  _±  1)  sln2 

N  +  1  [  2(N  +  3) 


The  first  integral  on  the  right  in  (3l)  is  evaluated  in  terms  of 
complete  gamma  functions,  and  the  second  integral  is  evaluated  by 
(30)  for  (j  =  (rr/2)  -  0,  "which  is  small  when  0  ~  rr/2.  Even  so,  if  N 
and  M  are  large  and  0  ~  tt/4,  it  does  no  good  to  carry  out  the 
interchange  of  N  and  M,  and  the  convergence  of  (30)  in  this  case 
would  be  slow.  Difficulty  (a}  would  be  met,  since  consecutive 
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terms  of  nearly  equal  magnitude  "but  opposite  sign  do  occur.  The 
associated  loss  in  accuracy  is  easily  seen  for  the  example: 

M  =  N  =  99,  0  =  tt/ 4 .  In  this  case,  1^2(50,50)  =  l/2,  the 

integral  on  the  left  side  of  (30)  is  equal  to  1.98  x  10"31.  The 
factor  sinN+1  @/(N  +1}  in  (30)  is  approximately  10-17,  and  thus 
the  second  factor  on  the  right  hand  side  of  (30)  must  he  of 

order  10” 14 .  But,  the  first  term  in  this  second  factor  is  unity, 
and  the  second  term  is  negative  and  greater  than  unity  in  absolute 
value.  Thus  the  second  factor  obviously  approaches  10” 14  neces¬ 
sarily  through  the  addition  of  nearly  equal  consecutive  terms  with 
opposite  sign. 

For  the  special  case  of  M  =  0,(b  =  l/2),  Fettis  sets  M  =  N 
in  (30),  and  uses  the  fact  that 

sinN  0  cosN  0  =  2_n  sinN  20 


to  derive 
0 

[  sinN  0  d0 


[2  sin  0/2]N+1 
N  +  1 


1  -  (N  -  1)  (N  +  1)  sln. 
2(N  +  3) 


2  (7; 


(e/2) 


(N  -  1)  (N  -  5)  (N  +  l) 

22  21  (N  +  5) 


sin' 


4  f  n 


(e/2) 


(32) 


where  0  =  20.  Equation  (32)  has  the  same  deficiency  for  large  N 
as  the  previous  relation,  (30).  The  values  of  0  range  from  zero 
to  tt/2,  and  again  nothing  is  gained  in  this  case  by  using  (31)  for 
0  near  tt/2.  If  M  =  0  in  (30),  then 


9 


sinN  0  d0  =  sinN+1  0  z 


(21)1 


sin21  0 


i=o  221  (il)2  (N  +  2i  +  1) 


.  (33) 


The  term  of  this  series  is  of  order  (/”3/2)  for  0  ~tt/2,  and 
again  (33)  would  not  lead  to  an  efficient  algorithm  for  such  0, 
even  though  the  terms  of  the  series  here  are  all  positive. 

M.  E.  Wise,  [17],  [18  ] ,  [19],  deals  with  the  inverse  problem, 
primarily,  of  finding  good  approximations  to  x,  the  percentage 
points  of  Ix,  given  Ix,  a,b.  Towards  this  objective,  he 
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advances  a  formula  for  Ix  in  [17]  which  is  derived  by  a  contour 
integration  in  the  complex  plane.  It  is  asymptotic  for  large  a 
and  b.  Wise  draws  attention  to  papers  by  E.  C.  Molina,  [8],  and 
C.  R.  Rao,  [ll],  and  another  of  his  own,  E 19 3 ,  in  which  similar 
results  are  derived  without  resorting  to  the  complex  plane.  Below 
we  give  our  own  derivation  of  Wise's  result.  The  equation  to  be 
derived  is  given  by  (43). 

Let 

t  =  e"u/N  ,Ost<co,OSxSl,  (34) 

where  N  ^  (a  +  b/2  -  l/2).  Then  substituting  (34)  into  (l)  and 
(2)  and  taking  their  difference  gives 


y-  -N  An  x 


0 


y 

=  i  2b~zJ^  e~U  [  sinh  (u/2W)]b”;L  du  t  (35) 

o 

The  term  [sinh  (u/2N)]b_1  is  expanded  in  powers  of  the  argument 
z  =  u/2N.  From  [1;  p  75,  equation  4.3.71] 

0° 

In  -A-ft-.?.  =  £  a^  z2n  ,  I  z  I  <  tt,  (36) 

z  n=i 

n  —  1  -p 

where  ,  and  B2n  is  the  2nrn  Bernoulli  number, 

n(2n)i 

[l;  p  810].  Hence 

sinb/b_1)  z  =  z13-1  exp  (b  -  l)  £  a^  z2n  .  (37) 

n=l 
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It  is  tedious  but  not  difficult  to  get  the  a^'s  for  small  values 
of  n.  Express  each  exp[(b  -  l)  a^  z2n]  in  its  power  series  about 
z  =  0  for  n  =  1,  2 ,  i  and  subsequently  carry  out  the  poly¬ 

nomial  multiplications.  The  first  six  a^  are  given  by 


a1  =  B2  =  i/6  , 


Thus 

exp[(b  -  1)8^  z2n]  =  Z  - —  fS-L  z2nl;  n  =  1,  2,  I  . 

i=o  il 

(39) 

The  product  over  n  is  taken  to  give 


exp  (b  -  l)  Z 
L  n=l 
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The  substitution  of  (40)  into  (35)  and  the  subsequent  term  by  term 
integration  of  (35),  with  respect  to  u,  gives 

Ba(a,b)  -  Bx(a,b)  a«  JL  P„(b)  +  T  (b  +  2) 

IT  L  24  •  N2 

+  -t5b  ~  7)  r  (b  +  4) 


27  •  32  •  5  *  N4  y 


1  zH  .(.35bj ±  95.1  r  (b  +  6)  +  ...1  , 

210  •  34  •  5  •  7  •  R6  y  J 


where 


ry(h)  =  /  e~u  u15'1  du  ,  h  >  0  , 


is  known  as  the  incomplete  gamma  function,  [1;  p  260]  .  The  final 
result  follows  by  dividing  both  sides  of  (42)  by  (2)j  thus  for 
4  =  3 

ix(a,b)  - 1  -  [ill  +  fo3  -  b)  ry(b  +  2A 

Nb  r(a)  Lr(b)  23  •  3N2  r(b  +  2) 

+  (b  -  1)  (5b  -  7)  Py(b  +  4) 

27  •  5  •  32  N4  r(b  +  4) 

+  1  -  1)  (55b2  ~  l^b  +  93)  +  6)  #  (43) 

(210  •  34  •  5  •  7)  N6  r(b  +  6)  -1 

The  ratio  ry(b  +  2i)/r(b  +  2i)  is  computed  by  the  recurrence 
relation 

r  (b  +  2)/r(b  +  2)  =  [r(b)/r(b)3 

-  e"y  (1  +  b  +  y)  yb/[b(b  +  l)  p(b)]  .  (44) 
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This  relation  is  derived  by  performing  two  integrations  by  parts  on 
Py(b  + 2 ),  as  defined  by  (42),  dividing  the  result  by  r(b  +  2)  and 
by  using  the  gamma  function  equation 

b  r(b)  =  r(b  +  1)  .  (45) 

The  series  given  by  (43)  is  quite  attractive  from  an 
asymptotic  standpoint.  Its  rate  of  convergence  from  some  numerical 
examples  appears  to  be  rapid  for  large  a  and  b. 

This  series  also  was  not  incorporated  into  our  program  for  a 
number  of  reasons.  First,  an  efficient  incomplete  gamma  function 
subroutine  is  needed.  Such  a  subroutine  does  not  seem  to  exist 
for  the  fast  computing  we  require.  Second,  the  storage  require¬ 
ments  for  such  a  routine  might  not  be  small.  Finally,  the  compu¬ 
tation  time  for  (43)  could  be  slowed  down  significantly  if  three  or 
more  terms  are  required,  because  of  the  cumbersome  nature  of  the 
coefficients.  Although  it  cannot  be  said  with  certainty,  since  the 
study  of  this  phase  was  very  limited,  it  appears  that  if  (43)  would 
be  more  efficient  than  the  method  we  employ  ( see  Sections  IV  and  V) , 
the  difference  would  not  be  impressive  as  far  as  computing  time. 
Certainly,  (43)  deserves  further  study. 

In  a  paper  by  I.  C.  Tang,  [14],  a  scheme  is  given  for  comput¬ 
ing  Bx(a,b).  The  basic  equation  in  his  paper  is  similar  to  (48),  (50) 
of  this  report.  He  has  developed  a  series  expansion  with  all  posi¬ 
tive  terms  in  place  of  the  usual  alternating  series  for  Bx.  The 
derivation  is  elegant,  however  the  relation  itself  was  known  to 
Soper,  [12].  Two  basic  problems  with  which  one  is  concerned,  for 
large  a  and  b,  in  the  application  of  Tang's  relations,  i.e.,  diffi¬ 
culty  { c ) ,  the  scaling  problem,  and  difficulty  [d],  the  computation 
of  starting  values,  are  not  discussed  in  [14]. 

This  section  is  closed  by  a  few  comments  on  two  digital  com¬ 
puter  programs  published  in  the  algorithm  section  of  the  Communi¬ 
cations  of  the  ACM. 
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W.  Gautschi,  [4],  describes  a  program  in  which  the  scaling 
difficulty  and  the  starting  value  problem  are  resolved.  The  basic 
relations  are  given  by 

Ix(a  +  n  +  l,b)  =  [1  +  (n+a+b  -  l)  x/(n  +  a)]  Ix(a  +  n,b) 

-  [(n  +  a  +  b  -  l)  x/(n  +  a)]  I  (a  +  n  -  l,b) 

(46) 

Ix(a,b  +  n  +  l)  =  [l  +  (n  +  a  +  b  -  l)  (l  -  x)/(n  +  b)]  Ix(a,b  +  n) 

-  [(n  +  a  +  b  -  l)  (l  -  x)/(n  +  b)]  Ix(a,b  +  n-l). 

(47) 

Nevertheless,  his  program  would  not  be  suitable  for  our  purpose 
when  a  and  b  are  large,  because  of  difficulty  {e).  For  example  if 
a  >  60  say  103  or  104  and  b  is  approximately  60,  it  would  require 
the  computation  and  summing  of  60  plus  10s  or  104  terms  of  (46) 
or  ( 47) . 

The  other  computer  algorithm  was  designed  by  0.  G.  Ludwig, 
[7].  It  was  programmed  at  NWL  for  the  ZBM  7030  (STRETCH)  by 
Mr.  Robert  Belsky  in  the  interim  period  of  development  of  the 
method  described  in  Section  IV.  Ludwig's  procedure  worked  quite 
well.  In  his  procedure,  four  sums  are  generated  in  every  case, 
whereas  in  the  present  method  no  more  than  four  occur,  and  in 
some  instances  only  one  summation  is  required,  e.g.,  when  b  is  an 
integer.  Moreover  if  x  >  l/2  and  a  is  large  Ludwig's  method  . 
requires  summing  over  approximately  the  integer  part  of  a  elements. 
This  leads  to  inefficiency  for  large  a  as  mentioned  previously. 

The  method  for  computing  Ix,  as  described  in  the  next  sec¬ 
tion,  was  developed  by  the  authors.  Although  it  includes  some 
relations  in  common  with  those  mentioned  in  some  of  the  preceding 
papers,  it  is  basically  a  complete  method  in  its  own  right,  since 
it  dispenses  with  all  the  difficulties  given  on  page  7  satis¬ 
factorily,  whereas  none  of  the  methods  described  in  this  section 
have  this  overall  feature . 
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IV.  AN  EFFICIENT  METHOD  FOR  COMPUTING  Ix(a,b) 


This  section  contains  the  main  results  of  this  report.  The 
analysis  that  was  developed  for  computing  Ix  is  separated  into 
3  cases  as  follows : 

A  -  a  or  b  is  a  positive  integer  no  greater  than  60 ; 

B  -  Neither  a  nor  b  is  an  integer,  and  a  ^  60; 

C  -  b  is  not  an  integer  and  a  >  60 . 

The  primary  ideas  or  motives  behind  the  method  are: 

(1)  that  a  and  b  can  be  represented  by  k  or  k-l/2  and  j  or 
j-l/2, respectively  where  j  and  k  are  positive  integers  such  that 
1  s  j  s  60,  1  £  k  ^  108; 

(2)  that  all  sums  will  be  finite  so  no  truncation  error 
occurs, with  two  exceptions;  in  these  cases  the  truncation  error  is 
rigorously  and  sharply  bounded  (See  discussion  on  (81)  and  the 
evaluation  of  An  r(s)  under  Section  V) j 

(3)  that  no  procedure  be  used  which  requires  summation  over 
k(a  =  k  or  k  -  l/2),  unless  a  £  60; 

(4)  that  no  alternating  power  series  are  evaluated. 

It  will  be  assumed  throughout  that  Ix(a,b)  is  to  be  computed 
to  an  accuracy  of  [ [log1Q  l/e]]  decimal  digits,  where  e  is  assigned 
and 


Case  A:  b 
If  b 


[[sj]  =  greatest  integer  in  s. 
j,  and/or  a  =  k  s  60  (See  Flow  Charts  (T)  ,(2)) 
j,  Ix  can  be  computed  from 


j 

Ix(a^)  =  2  a*  , 


(48) 
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where 


a±  =  x a  r(a  +  (1  -  x)1"1  ,  (49a) 

T(a)  r(i) 


ai  =  xa  =  Ix(a,l)  . 

If  a  =  k  *  60,  I*  may  be  computed  by  using  (5),  i.e.. 


(49b) 


I  (a,b)  =  1  -  Ix(a,b)  =  1  -  Z 

1=1 


(50) 


where 


and 


bj.  =  (1  -  *)b 


r(b  +  i  -  1)  i-1 

r(b)  r(i) 


Ix(a,b)  s  I1_x(b,a)  . 


(51) 


(52) 


The  choice  between  (48)  and  (50)  is  made  accordingly: 

if  a  ^  k  (a  not  an  integer),  b  =  j,  then  use  (48); 

if  a  =  k  2  60,  b^j  (b  not  an  integer),  then  use  (50); 

if  a  =  k  s  60,  b  =  j,  then  ■use  (48)  if  j  s  k  and  use 

(50)  if  j  >  k. 

The  derivation  of  (48)  is  given  in  Appendix  A.  Equations  (50), 
(51)  follow  directly  from  (5)  and  (48). 


The  remainder  of  the  analysis  on  Case  A  will  be  with  respect 
to  (48)  since  the  results  for  (50)  are  analogous  by  the  substitu¬ 
tions  implied  by  (5). 

A  complication  arises  from  (49a)  because  of  the  gamma  func¬ 
tions.  Although  each  a1  .must  remain  less  than  1  -  e  (otherwise, 
since  all  a.]_  >  o,  1,^1-  e)  the  individual  quantities  r(a  +  i  -  1) 
and  r(a),  and  even  their  ratio  r(a  +  i  -  l)/r(a),  can  exceed  the 
value  of  the  largest  single  precision  number  the  computer  can 
operate  on.  The  same  problem  is  manifest  in  the  b1  and  coef¬ 
ficients  given  by  (51)  and  (62),  respectively.  This  difficulty  with 
the  a1(and  bq)  is  resolved  by  the  following  scaling  procedure: 
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Let 


an  =  max  a.±  , 


then 


.min  j[[(a  -  l)(l  -  x)/x]]  +  1, 

n  = 

1 


The  result  given  by  (53)  is  easily  deduced  since,  by  the  defini¬ 
tion  of  an,  it  is  required  that 


an-i  <  an  )  2  ^  n  ^ 

(54a) 

an+i  ~  1  ^  n  ^ 

J  -  i. 

(54b) 

Inequalities  (54)  imply 

n  ^  [(a  -  l)(l  -  x)/x] 

+ 1  , 

(55a) 

n  ^  [(a  -  l)(l  -  x)/x ] 

? 

(55b) 

from  which  (53)  follows.  Inequalities  (55)  also  imply  that 
there  are  at  most  two  an  and  if  so  they  are  consecutive . 


Having  found  an  expression  for  an,the  in  an  can  be  computed 
by 


k/l 
k  =  1  . 


(53) 


in  an  =  a  in  x  +  (n  -  1)  in(l  -  x) 

+  in  r(a  +  n  -  l)  -  in  r(a)  -  in  r(n) 


(56) 


Various  sensings  are  made  on  £n  an  from  which  it  may  be  usually 
concluded  if  Ix  £  e  or  Ix  ^  1  -  e.  Thus  all  the  are  under 
control,  at  this  stage,  since  none  is  larger  than  an  =  exp  [jgn  an] 
which  must  remain  less  than  1  -  e  as  explained  above.  The  pro¬ 
cedure  is  brought  forth  in  detail  in  Flow  Chart  (2,  . 
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The  ai(i  ^  n)  are  computed  by  the  following  extremely  simple 
and  efficient  recurrence  relations: 

ai+i  =  (—-"!■■--)  (1  -  x)  a.  ,  nsisj-l,  (57) 

ai  =  (t - - - r)  f  -  1  ]  ai+i  >  1  <:  i  <  n  -  1  ,  (58) 

\i  +  a  -  1/  \1  -  xj 

which  are  easily  derived  from  ( 49a) . 

The  computation  of  (56)  requires  a  method  for  evaluating 
in  r(s)  directly,  where  s  is  used  to  represent  the  argument  of  any 
natural  logarithm  that  appears  in  this  report.  The  method  by  which 
this  is  accomplished  is  described  in  Section  V,  and  by  flow  charts 

(?)  and  Q). 

Case  A  is  concluded  by  noting  the  following  advantages: 

(1)  All  terms  of  the  sums  in  (48)  and  (50)  are  of  like  sign. 

(2)  The  series  to  be  summed  are  finite  series  with  the  number 
of  terms  to  be  summed  not  exceeding  60.  Thus  no  truncation  error 
need  occur  (actually  one  is  introduced  by  a  sensing  in  the  program 
which  permits  truncation  of  the  series  if  any  of  the  a.,  or  b± 
become  less  than  specified  tolerances.  See  Flow  Chart  (?).). 

(3)  The  magnitudes  of  the  aj  are  kept  under  control  for 
any  k  such  that  l/2  £  a  ^  10s. 

(4)  The  procedure  is  efficient. 

Case  B:  a  s  60,  a  =  k  -  l/2,  b  =  j  -  l/2.  See  Flow  Charts  (1^, 

©  >  (i-  • 
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In  this  case  Ix(a,b)  is  computed  from 


Ix(a,b)  =  Ix(a,l/2)  +  •/l~rx'  E  xa  Aa,+  j  -  ll^L  (l  -  x)1'1 

X  i=1  r(a)  r(i  +  1/2)  (59) 


ix(a,i/2)  =  ix(i,  i)  -  m  sr^n  kz 


i=1  rU  +  1/2)  r(i/2] 


x1"1  (60a) 


I v(—,  b')  =  Ix(i,  +  /x  /I  -  x  Z 


L=l  r(i  +  1/2)  rC  1/2) 


(1  -  x)1-1 


(60b) 


such  that  the  arc  tangent  is  between  0  andTT/2.  The  derivations 
of  (59),  (60)  are  given  in  Appendix  A. 


We  introduce  the  notation 


;  =  xa  F(a  +  t  -  U^L  (l  -  x)1-1/2 

r(a)  r(i  +  1/2) 


di  r(i  +  1/2)  r(i/2) 


such  that  (59)  becomes 


Ix(a,b)  =  —  tan' 
TT 


—  tan“ 1  (  x [  -  ,/x" 
tt  V/l  -  x  / 


_  k-i  j-l 

x  E  d.  +  E  c, 
1  i  *1 
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The  terms  di  are  generated  by  the  simple  recurrence  relation 


di«  * x  \~i I  ai  -  i s  1  s  k '  2  ’ 


where 

d,  =  “  • 

1  TT 

Wo  scaling  problem  occurs  with  the  d±  terms* 

The  c±  terms  require  scaling;  it  is  done  in  the  same  way 
that  the  a ^  were  scaled. 


c n  =  max  c^  , 


then 


n  =  min  i[[(a  -  l)(l  -  x)/x  +  -]]  }  j  -  ll  (66) 


where  the  result  is  derived  from 


*n+i  _  n  +  a  -  1/2 
Cn  n  +  1/2 


(l  -  x)  £  1  , 


n  +  a  -  5/2 
n  -  1/2 


(1  -  x)  ^  1  ; 


the  double  square  bracket  notation  was  defined  on  page  19.  It  is 
known  that  cn  is  less  than  1  -  Ix(a,l/2)  as  otherwise  Ix(a,b)  is 
equal  to  one.  From  this  point,  the  scaling  proceeds  exactly  as 
for  the  an.  The  recurrence  relations  for  the  c.  are: 

ci+1  =  (1  "  x)  Ci  >  n  *  i  £  j  -  2  ,  (69) 


i  +  1/2 
i  +  a  -  l/2 


i+1  > 


1  £  i  £  n  -  1  .  (70) 
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The  ci  are  also  summed  in  the  same  order  that  the  a1  were,  i.e.,  the 
ci  with  i  >  n  are  computed  and  summed  in  increasing  order  of  i.  Then 
the  remaining  c±  are  computed  and  summed  in  decreasing  order  beginning 

at  i  =  n  -  1.  Refinements  in  the  program  which  are  not  included  here 
may  be  gleaned  from  Flow  Charts  (3)  and  (4)  *  The  favorable  factors 
listed  for  Case  A  on  page  22  also  apply  for  Case  B,  with  the  excep¬ 
tion  that  two  rather  than  only  one  sum,  with  as  many  as  sixty  terms, 
may  have  to  be  evaluated. 

Case  C:  a  >  60,  b  =  j  -  l/2.  See  Flow  Charts  (T)  ,  (3)  ,  (5)  * 

Case  C  is  by  far  the  most  difficult  of  the  three  cases  to  evalu¬ 
ate  Ix(a,b).  The  beta  ratio  is  again  given  by  (59),  however  (60) 

cannot  be  used  for  the  computation  of  Ix(a,l/2)  because  the  summa¬ 
tion  in  (60a)  runs  over  k,  where  this  integer  can  be  much  larger  than 
60.  Thus  the  problem  here  reduces  to  finding  an  efficient  procedure 
for  evaluating  Ix(a,l/2)  when  a  is  large.  After  considering  some  of 

the  methods  proposed  in  the  literature,  [3],  [12],  [ 18 ] ,  [20],  it  was 
decided  to  proceed  by  an  entirely  different  approach,  that  of  using 
Gaussian  quadrature,  [6;  p  319].  This  technique  was  chosen  because 
the  truncation  error  E  (^  e ')  could  be  sharply  and  rigorously  bounded, 
and  moreover  the  error  bound  indicated  that  a  surprisingly  low  order 
Gaussian  formula  would  suffice  for  the  accuracy  desired.  The  details 
of  the  critical  steps  in  the  proofs  required  for  the  bound  E'  of  the 
Gaussian  error  term  are  relegated  to  Appendix  B,  otherwise  the 
analysis  needed  follows.  We  begin  with  some  preliminaries. 

Apply  the  transformation  t  =  1  -  u2  to  Bx(a,l/2),  so  that 

1 

Bx( a, l/2)  =  Z  f  (1  -  u2)a-1  du  .  (71) 


The  following  notation  is  introduced: 
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One  can  then  write 

Ix(a,l/2)  =  lx(a,l/2;X)  +  I1_x2(a,l/2;1) , 


(74) 


with  the  objective  of  making  the  last  term  in  (74)  small,,  i.e.,  for 
a  given  €/;  >0  to  rind  a  X  such  that 

I1_x£(aU/2J1)  £  e".  (75) 

A  function  X(e")  which  satisfies  (75)  is  given  by  (80).  The  derivation 
follows.  From  (73)  . 
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where  z  =  /"a  -  1  u.  But  the  last  integral  in  (77)  is  negligible  for 
a  ^  60, since 


J  e"u2  du  .  (See  [l],  p.  298). 

v/a-l 


Hence 


^^(aA/s;!) 


1) 


a 


1) 


(- 


{£.  _ -X2(  a-l) 

/a  -  1 


=  - ^  ■  exp  (a  -  l)  (-  -z  X2l/i)  1 

M  /a-l  L  V  1  '/J 


(i  -  x2)^1  * «# 

M  /a  -  1 


(78) 


By  solving  (78)  for  X  one  obtains  that 


X  ^ 


(79) 


Inequality  (79)  is  relatively  sharp;  it  can  be  slightly  improved  provided 
one  is  willing  to  solve  a  transcendental  equation  for  X  and  accept  the 
corresponding  increase  in  computing  time  which  would  result. 
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The  smallest  value  of  X  is  chosen  from  (79)  for  x(0  so  that 


X(e")  = 


r  /  ^ .  a-in 

i/2  I 

-1 

1  /M  /a  -  1  c//\ 

[-  rh e  J 

L  '  yw 

1/2 


(a  -*  “) .  (80) 


If  a  is  large  x(e")  =  X  should  he  determined  from  its  asymptotic 
fom  as  given  in  (80).  One  can  now  deduce  from  (78)  and  (80)  that 
the  upper  limit  (unity)  of  the  integral  in  (71)  can  be  replaced  by 


the'  smaller  quantity  x(e")  =  X. 
value  of  Ix( a,  l/2)  is  less  than 

from  (73)  and  (74). 


If  X  is  less  than  /I  -  x  ,  then  the 
e" ;  a  fact  that  is  easily  concluded 


Having  dispensed  with  these  introductory  results,  the  basic 
objective  here  of  deriving  a  truncation  error  bound  for  the 
Gaussian  integration  procedure  is  now  carried  out. 


The  exact  error  term  as  a  result  of  using  Gaussian  quadra¬ 
ture  of  order  m,  0(m),  [6;  p  324],  to  numerically  compute  the 
integral  of  a  function  f(t),  with  a  sufficient  number  of  denva 
tives,  over  [-  1,  l]  is  given  by 


22m+l  (m.)4  )  -  1  <  tx  <  1  ,  (81) 

(2m  +  1)  [(2m)!]3 


where  f(2m)(t)  means  the  2mth  derivative  of  f(t)  with  respect  to 
t.  The  integral  of  (73)  is  transformed  so  that  the  limits  of 
integration  become  -  1  and  1.  The  usual  transformation 


u 


(x  -  /r^~x)  t  +  (x  +  /!  -  *1 

2  2 


(82) 


applied  to  (73)  gives 


Ix(a,l/2;X) 


X  -  /I  -  X 
M 


F(u)dt  . 


-  1 


(83) 
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Similarly  by  applying  the  transformation 

u  =  (v/1  -  x/2)  (l  +t)  (84) 

to  Ix(a, l/2),  defined  in  (52),  the  result  is 

l 

Ix(a,l/2)  =  &  - . J  F(u)dt  ,  (85) 

-l 

where  F(u)  represents  the  Integrand  in  (73).  It  is  important  to 
consider  I  here  as  well  as  I  .  The  total  interval  of  integration 
as  specified  in  (73)  is  (x  -  /I  -  x),  however  if  we  apply  (83)  only 
when  x/2  <  /I  -  x  (and  (85)  only  when  x/2  s  /I  -  x) ,then  essentially 
the  total  integration  interval  is  never  larger  than  x/2  or  half  the 
.maximum  value  of  (X  -  /I  -  x) .  This  leads  to  a  decrease  in  E '  by  a 
factor  of  2+(2m+1),  since  the  integration  interval  appears  in 
f(£m)(t)  explicitly  to  the  (2m+l)  power  (see  Equation  87).  The  term 
x/2  is  obviously  never  larger  than  l/2  and  generally  will  be  quite 
small.  For  example,  if  a  =  104,  e"  =  9  x  10 then  x/2  =  .024 
from  ( 80 ) . 

The  2mth  derivative  of  f  with  respect  to  t  is  needed.  The 
integrand  from  (83)  and  (85)  is  given  by 


(86b) 
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Therefore,  indicating  the  2mth  derivative  with  respect  to  u  of 
F(u)  by  F^2m)  ,  f  ^2ra)  (ti)  in  (81 )  is  given  by 


2  / X  -  /I  -  x 
M\  2 


2  / /I  «•  x 
M  \  2 


FS2m)(Ul)  , 


F^UJ  , 


\  <  '/l  -  x  , 


(87a) 


(87b) 


where  u  =  corresponds  to  t  =  ti,and  it  is  understood  tj_  (and 
also  u-J  is  different  in  (87a)  from  (87b).  The  effect  of 

reducing  the  integration  interval  is  evident  in  (87).  It  is 
observed  that  the  second  factors  on  the  right  hand  side  of  (8?) 
are  bounded  by  (x/4)2m+1.  The  principal  result  we  wish  to  derive 
is  the  following  expression  for  E '  , 

E  £  E'  3  §  r-UZgJg-mtlC?-OLl  [ - r(a).J  e/  (88.) 

_(2m  +  l)  [ ( 2m) 1 3 2 J  p( a  -  m) 

subject  to  the  constraint  that  a  -  1  >  2m  +  l/2.  Since  a  >  60 
here  and  m  will  turn  out  to  be  in  the  neighborhood  of  ten,  the 
constraint  is  insignificant  for  our  application. 


u  3  [(1  -  u2)n)  ,  o  £  u  *  1  , 


so  that 


U  =  F  ^2m* 

a-i  2m  u 


It  is  shown  in  Appendix  B  that 

[[r/a]] 

U„  -  z  (-  p*-1  — .  -gr"21 


rl r(n  +  1) 


,2\n-r+I 


il(r  -  2i)l  p(n  -  r  +  i  +  1) 
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and  also  that  Un  r  satisfies  the  following  ordinary  differential 
equation 


(1  “  u2)  U^r  +  2(n  -  r  -  l)u  p  +  (r  +  l)(2n  -  r)  p  =0 

(92) 


for  positive  integers  r  and  real  numbers  n  ^  r. 

The  key  idea  which  leads  to  a  useful  bound  on  Un;  r  is  that 
the  absolute  values  of  the  extrema  of  Un^  r  form  a  decreasing 
finite  sequence  on  [0,1]  for  n  >  r  +  l/2.  Two  closely  related 
proofs  of  this  statement  are  given  in  Appendix  B.  It  is  the 
crucial  step  in  the  sequence  of  steps  employed  for  bounding  E. 
Thus  assuming  the  statement  true,  it  follows 

|un,  r(u)l  lUn,  p(°)l  f  r  even-  (93) 

For  n  =  a  -  1,  r  =  2m,  one  obtains  from  (91) 


lu  (0)1  =  IgsliEL&L 

a-h2®'  m I  r( a  -  m) 


(94) 


and  the  desired  result  given  by  (88)  follows. 

The  graphs  on  pages  36,  37  contain  curves  of  (-log1D  E7) 
versus  a  based  on  (80)  and  (88),  for  given  values  of 
e"  and  m.  Their  purpose  is  to  indicate  the  smallest  order  of 
Gaussian  integration,  0(m),  which  can  be  used  for  a  given  e', 
where  e'  represents  the  upper  bound  on  E'.  The  results  as 
graphically  set  forth  clearly  substantiate  the  remark  .made 
earlier  that  very  low  order  Gaussian  integration  formulas  will 
suffice  for  the  evaluation  of  I  (a,  l/2)  for  large  a.  For 
example,  in  the  computing  program,  as  it  is  now  operating  m  =  10 

is  used  with  e"  =  9  x  10 _11,  e'  =4.5  x  10 _11,  and  it  is 
apparent  from  the  graphs  that  this  value  of  ,m  is  adequate  for 
all  a  e  [60,10s] . 
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The  procedure  by  which  the  curves  were  constructed  for  a 
given  e"  was  as  follows: 

(a)  a  sequence  of  positive  integers  was  chosen  to  represent 
various  0(m) , 

(p)  xCe")  was  then  computed  by  (80)  for  a  sequence  of  values 
a  e  [60,i04°],  and  a  given  value  of  e", 

(7)  these  computed  values  of  X  with  their  corresponding  a  values 
were  then  used  in  (88)  to  compute  log1  E'. 

Thus  the  0(m)  to  be  used  for  a  given  e  '  can  usually  be  estimated 
conservatively  from  the  graphs.  A  precise  0(m)  can  always  be  deter¬ 
mined  by  computing  a  set  of  X  from  (80)  and  the  associated  E'  from 
(88)  for  various  m  and  a.  One  observes  that  generally  logl0  E '  is  a 

very  slowly  increasing  fuction  of  a. 

This  section  is  concluded  with  the  explicit  formulas  used  for 
the  Gaussian  quadrature  of  Ix(a,l/2).  They  are: 


Ix(a,l/2)  -  (X 


m 

♦  Z 

1  =  1 


(X  -  ./I  -  x) 


1  + 


+ 


yr 


X 


a-  1 


J 


X2  <  4  (l  -  x)  , 


Ix(a,l/2)  fev/T 


r(a  +  1/2) 

r(a)  r(i/2) 


m 

•  Z 

i  =  l 


Wi  /  1 


X2  ^  4  (l  -  x)  , 


(95) 


(96) 
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where  the  y1  and  wi  are  the  Gaussian  abscissae  and  weights  , 
respectively,  of  order  m,  0(.m),  [l;  p  916].  Since  the  last  term  in 
(74)  is  always  non- negative  and  no  larger  than  e"  for  \  which 
satisfies  (80),  it  follows  that 

Ix(a,l/2)  -  [lx(a,l/2;X)  +  e"/2] 

since 

Ix(a,l/2)  ^  Ix(a,l/2;X)  . 

This  accounts  for  the  additional  e"/2  in  (95). 

The  Wj_  and  (l  +  y.  )/2  are  tabulated  below  for  0(10)  to  14 
significant  digits  on  f-1,  1],  where  y^  are  the  Gaussian  abscissae 
and  the  Gaussian  weights . 


2  e'/2  ,  (97) 

(98) 


(1  +  yi)/2 

0.0130  4673  5791  414 
0.0674  6831  6655  507 
0.1602  9521  5850  49 
0.2833  0230  2935  38 
0.4255  6283  0509  18 
0.5744  3716  9490  81 
0.7166  9769  7064  62 
0.8397  0478  4149  51 
0.9325  3168  3344  49 
0.9869  5326  4208  59 


0.0666  7134  4308  688 
0.1494  5134  9150  58 
0.2190  8636  2515  98 
0.2692  6671  9310  00 
0.2955  2422  4714  75 
0.2955  2422  4714  75 
0.2692  6671  9310  00 
0.2190  8636  2515  98 
0.1494  5134  9150  58 
0.0666  7134  4308  688 


The  flow  Chart  (5)  covers  this  part  of  the  program. 
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The  use  of  (59),  with  Ix(a,l/2)  precomputed,  as  above,  gives 
Ix(a,b). 


The  quantities  e,  e' ,  and  e"  which  appeared  in  this  section 
are  briefly  ^ummarized.  The  number  e  is  specified  slightly  less 

than  5  x  10  ^+1  where  p  is  the  number  of  decimal  digits  to  which 

Ix(a,b)  is  to  be  computed.  The  number  e/  is  specified  and  is  used 
for  bounding  the  truncation  error  due  to  Gaussian  quadrature  which 
is  used  to  evaluate  Ix(a,l/2)  when  a  ^  60.  The  graphs  on  pages 
36-37  are  a  guide  to  determine  the  0(.m)  for  given  e ',  e". 

Generally  e'  is  taken  equal  to  e/4.  The  number  e"  is  taken  equal 
to  2(e  -  e')/3.  This  number  is  used  in  (59)  to  bound  the  cn.  The 
details  are  shown  in  Flow  Chart  (§)  .  The  quantity  e"  is  also  used 
in  (80)  to  reduce  the  Gaussian  interval  of  integration  from 

[/l  -  x,  l]  to  [/l  -  x,  X]. 


The  e- quantities  are  used  primarily  to  insure  that 
Ix(a,b)  is  computed  within  e  when  (59)  is  used.  Thus  if 

e '  -  e/4  and  e"  =  2(e  -  e')/3  =  e/2,  the  following  analysis 
shows  that  the  required  accuracy  is  attained. 

Ix(a,b)  =  Ix(a,l/2jX)  +  I1_xa(a,l/2jl)  +  S  ,  (99) 

where  S  denotes  the  summation  of  j-1  terms  in  (59);  and  a  similar 
relation  holds  between  the  computed  quantities,  distinguished  by 
asterisks  from  the  corresponding  true  values  in  (99).  Then,  taking 
differences  and  using  the  triangle  inequality, 

|lx(a,b)  -  l£(a,b)|  §  |lx(a,l/2;X)  -  I*(a,l/2;X)| 

+  | I1_^2(a,l/2;1)  -  I^_x2(a,l/2;1)|  +  A  -  x  |s  -  S*|  , 

(100) 

where  the  value  given  to  I^_x2(a,l/2;  l)  is  explained  below. 
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But  the  first  term  on  the  right  in  (100)  does  not  exceed  e7  9 
or  e/4,  through  the  choice  of  the  proper  order  of  Gaussian 
integration,  how  0  =  I1^2(a^l/2 ;1)  -  e" ,  by  (78).  Hence , 

reasoning  as  in  (97)  and  (98),  we  arbitrarily  take  1^  -xa(a,l/2jl) 
as  e." /2.  =  e/4,  and  this  guarantees  that  |  I1_-^2(a,l/2;l) 

-  I*_^2(a,l/2 jl) |  -  e/4.  The  last  term  in  (100)  does  not  exceed 

e"  or  e/2,  as  shown  by  the  .method  of  determining  the  number  of  terms 
computed  in  the  summation  (see  Flow  Charts  (3)  and  (4)).  Thus 
|l  (a,b)  -  I*(a,b) |  =  e/4  +  e/4  +  e/2  =  e,  as  was  to  be  shown. 


The  program  is  presently  set  for  obtaining  Ix(a,b)  to  within 

two  units  in  the  tenth  decimal  digit  for  a  ^  10s.  The  e-quantities 
are  specified  by 

e  =  1.8  X  10 -1°,  e '  =4.5  x  10 _11,  e"  =  9.0  x  10 "1X. 
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V.  COMPUTATION  OF  to  r(s),  K  s  In  r(a  +  c)  -  in  r(a) 

The  cases  A,  B,  C  which  have  been  described  above  require  the 
computation  of  in  r(s)  or  K  to  high  accuracy,  where  s  represents  a 
positive  integral  multiple  of  one-half  and  c  takes  the  values  n  -  1 
or  l/2.  This  is  dealt  with  in  an  efficient  manner,  at  the  expense 
of  two  hundred  storage  locations,  by  storing  the  value  of  in  r(s) 
for  s  =  l/2(l/2)l00  to  the  full  accuracy  of  a  single  precision  num¬ 
ber  (which  is  fourteen  significant  digits  on  STRETCH)  and  by  using 
asymptotic  series  for  in  r(s)  or  K  when  s  >  100. 

In  such  cases,  it  would  seem  convenient  to  always  use  the 
asymptotic  series,  for  in  r(s),  [l;  p  257],  which  is  given  by 


in  r(s)  -  (s  -  l/2)  in(  s  -  l)  -  (s  -  l)  +  (l/2)  in  2 rt 


1  1  1  1  1  1 
+  12  s  -  1  "  360  (s  -  l)3  +  1260  ( s  -  l)5 


.  (101) 


where  the  sum  of  the  first  five  terms  is  sufficient  for  thirteen  decimal 
digit  accuracy  with  s  >  100.  It  is  observed  however  that  in  every  case 
where  in  r( s)  is  needed  actually  the  difference  K  appears .  The  use  of 
(101)  to  compute  the  two  logarithmic  terms  of  K  separately  leads  to  a 
prohibitive  loss  of  significant  digits  if  a  is  very  large.  This  may 
be  seen  by  observing  that  the  dominant  term  in  (101)  for  s  =  a  +  c  or 
a  is  of  the  order  of  a  in  a.  Thus,  upon  subtraction,  the  undesirable 
loss  of  digits  occurs,  e.g.,  if  a  =  104  and  c  =  l/2  four  digits  are 
lost.  If  a  =  10s,  c  =  1/2  then  in  r(l08  +  l/2)  -  in  r(l08) 

=  1742068075.3142  -  1742068066.1038  =  9.2104,  so  that  in  this  case  nine 
digits  are  lost. 

This  difficulty  is  resolved  by  introducing  the  following  asymptotic 
series  for  K,  if  a  >  100, 


in  r(a  +  c) 


in  r(  a)  as  c  -  —  — 
2  a 


in(l  +  c/a) 

— 


1 


1  c 

2  a 


+ 


c 


in(  a  +  c)  - 


12 


1 

a 


1 

a  +  c 


‘1 

1  1 

1 

'l 

1 

_a3 

(a  +  c)s_ 

1260 

(a  +  c)8_ 

(102) 
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The  series  can  be  derived  by  the  use  of  the  standard  Stirling  approxi¬ 
mation  to  in  r(s).  The  first  expression  in  square  brackets  on  the 
right  hand  side  of  (102)  is  evaluated  by  the  series 


fn(l  +  d}  _ 

e 


y  + 


2vg 
2  +  0 


(103) 


where  0  =  c/a ,  and  y  =  0/(2  +  @) .  The  series  in  (103)  is  easily  and 
efficiently  generated  by  the  recurrence  relation 


^ "  H§  1 3) y2  Ap  - 1  *  P  ~  ^  2> 


(104) 


where 


Ap  h  y2P-  (2p  +  3)  ,  A0=l/3. 


(105) 


Either  five  or  ten  terms  of  this  series  is  used  to  attain  fourteen 
digit  accuracy  such  that  if 

1.  0  <  9  ^  0.15  five  terms  are  used; 

2.  0.15  <  0  <  0.6  ten  terms  are  used. 

It  is  also  necessary  to  retain  the  series  given  by  (101)  for  those 

cases  when  the  value  of  in  r(a  +  c)  is  not  stored  and  yet  a  <  100,  e.g., 

if  a  +  c  =  140  and  a  =  90,  In  such  cases  no  significant  loss  of 

digits  will  occur  in  computing  the  two  logarithmic  terms  of  K 

separately. 

Thus  if  a  +  c  ^  100 ,  K  is  obtained  by  table  look-up.  If 
a  +  c  >  100  but  a  <  100  then  in  r(a  +  c)  is  computed  from  (101)  and 
in  r( a)  by  table  look-up.  If  a  >  100  then  K  is  computed  from  (102) 
and  (103).  The  details  are  given  in  flow  charts  (6^ \9  (?)  . 
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VI.  COMPUTING  PROGRAM  FOR  Ix(a,b)  -  FLOW  CHARTS 


The  numerical  calculation  of  Ix(a,b)  for 

1  ^  a  ^  10s  .  I  £  b  ^  60  , 

2  2 

a  =  k,  or  a  =  k  -  ±  ,  b  =  J,  or  b  =  J  |  , 

where  j  and  k  are  positive  integers,  is  based  on  equations  (48), 
(50),  (64),  (80),  (95),  (96),  (97),  (lOl),  (102),  (103)  of  the 
last  two  sections.  The  program,  as  outlined  on  the  Flow 
Charts  (l)  -  © ,  has  been  coded,  as  a  subroutine,  for  the  NORC  and 
the  IBM  7030  (STRETCH)  in  absolute  machine  language,  Mr.  Travis 
Herring  prepared  the  STRETCH  coding.  . 

The  inputs  to  the  program  are  a,  b,  x,  e,  e'.  The  e- quanti¬ 
ties  are  discussed  on  pages  34-55  of  the  last . section.  If  the  number 
of  decimal  digits  required  in  Ix(a,b)  is  other  than  2  units  in  the 
tenth  decimal  place,  then  this  could  necessitate  a  change  in  the 
number  of  Gaussian  .multipliers  required  as  a  result  of  a  change 
in  e'  and  e". 

Two  constants  appear  in  the  flow  charts  which  depend  on  €. 

They  are  identified  by  the  letters  f  and  g  and  are  defined  by 


f  =  in  e"  =  in  2(c  -  e  ')/3 
g  =  in  € 


(106) 


Generally,  the  notation  used  in  the  flow  charts  allots  lower 
case  letters  to  numerical  values  and  identifies  the  machine  location 
in  which  that  value  is  stored  by  the  corresponding  upper  case  letter; 
thus  the  storage  location  for  the  number  o  would  be  Z. 

There  are  a  total  of  seven  flow  charts  starting  with  the 
master  flow  chart  in  which  the  over-all  computing  procedure  is 
outlined. 

The  average  computing  time  on  STRETCH  is  2.6  x  10“3  seconds 
per  case.  This  would  mean  that  an  average  computing  time  per  case 
on  an  IBM  7090  would  be  about  8  milliseconds.  The  average  time  on 
STRETCH  was  determined  by  running  large  sets  of  cases  for  random 
choices  of  x  values.  Also  very  many  cases  were  run  by  taking  equal 
increments  in  the  variables  a,  b,  x  and  essentially  spanning  the 
space  generated  by  these  variables. 
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It  is  easily  observed  from  Flow  Charts  ©  and  ©  that  if 
values  of  b  higher  than  sixty  are  desired  then  the  number  of  terms 
to  be  summed  in  such  equations  as  (48)  and  (64)  are  correspondingly 
increased.  If  b  were  made  excessively  large  the  procedure  given 
here  would  be  inefficient . 


For  easy  reference,  the  basic  formulas  used  in  the  program 
are  given  again  here  with  the  same  equation  number  they  carried 
previously. 

Ix(a,b)  =  1  -  I,__(b,a)  =  1  -  Ix(a,b)  ,  (5) 


Ix(a,b)  =  E  xa  ~  ~~  *■  (1  -  x)  1  ,  b  =  J  .  (48) 
l=l  r(a)  f(i) 


IT(a,b)  =  I  (a, 1/2)  +  E 


E  xa  ■-  j-  - ■  1//g,\  (1  -  x)1"1/2  , 

1=i  r(a)  r(i  +  1/2) 


a  =  k  -  —  i 
2  ! 


ib  =  *"2 


I  (a,l/2)  =  l  (l/2 ,  1/2)  -  Sk  /T^~k  E 


=1  r(i  +  1/2)  r(i/2) 


b  =  l/2,  a  =  k  -  l/2  £  60 


(60a) 


Ix(l/2,b)  =  Ix(l/2,  l/2)  +  Sk  /l  -  x  E 


1=1  r(i  +  1/2)  r(i/2) 


(l-x): 


a=i,b=J-i.  (60b) 
2  2 
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xK)  =  [- 

Ix(a,l/2) 

Ix(a,l/2) 

in  r(  s)  5; 

+ 

in  r(  a  +  c 


(61) 


M  ,/a~-  1 

/rT 


°°).  (80) 


(X  -  /"l  -  x)  +  ,rl/^ 

r(a)  r(i/2) 


X2  <  4  (1  -  x)  . 


a  /T~^~x  Ila-t  1/2), 
r(a)  r(i/2) 


r 

1 

{1- 

2| 

V  2  / 

[ 

J 

(95) 


X2  £  4(1  -  x)  . 


(96) 


(s  -  l/2)  in{  s  -  1)  -  ( s  -  l)  +  (l/2)  in  2tt 


1  1 


1  1  +  1  1 


(101) 


12  s  -  1  360  (s  -  l)3  1260  (s  -  l)' 


•  •  • 


:)  -  in  r(a)  2  c--- 

2  a 


ln(l  +  c/a)  , 

.  c/a  ‘  . 


—  —  +  c  in(a  +  c) - i 

2  a  12 


La  a  +  c 


(102) 


1 

1  1  1 

.1 

1 

H 

1 

h| 

1 

360 

-a3  (a  +  c)3  . 

1260 

La5  (a  +  c)5  J 
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ln(l  +  9) 
0 


1  =  -  y  + 


2yc 


2  + 


i  +  z!  +  z!  +  z!  + 

3  5  7  9  + 


=  ^2p~+~%J  y2  %>  -  1  >  p  - 1>  2)  *  ' 

where 

Ap  =  y2P/(2p  +3)  ,  A0  =  1/3  . 


,  (103) 


(104) 


(105) 
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MASTER  FLOW  CHART  Xx(a,  b) 


Ix(o,  b)  =  l-I|_x(b,  a)  s  l-Ix(a,  b) 

k  AND  j:  POSITIVE  INTEGERS 


INPUT:  a,  b,  X 
e.  «' 


1/2  <  b  <  60 


1/2  <  a  <  I08 


0  <  X  <  I 


b  =  j  OR  j  -  1/2 


a  =  k  OR  k  -  1/2 


OUTPUT:  Ix(a,b)  0<Ix(a,  b)<l 

PERTINENT  EQUATIONS  ARE  GIVEN  ON  PAGES  41  -  43 


44 


45 


START 


COMPUTATION  OF  EQUATION  (59) 
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m 


COMPUTATION  OF  EQUATION  60A 


Ix(a,  1/2)  =  Ix(l/2,  1/2)  /X  v/i-X  2  p(j+|/2)p[|/2)  X  »  Tx(l/2,  1/2) 


COMPUTATION  OF  EQUATIONS  (95),  (96) 


START 


EXIT  TO  EO.  59 
IF  b  ±  1/2 
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© 

COMPUTATION  OF  [loger(s)],  EQUATION  101 


*  VALUES  OF  loge  T(s)  STORED  IN  CONSECUTIVE  LOCATIONS  U|,  U2, -.UgOO 
ORDER  OF  s  FOR  s  =  1/20/2)100,  e  g.,  U|40  CONTAINS  log  TITO). 

1/2  loge  2 IT  =  0.9189  3853  3204  87 

n3  IS  NOT  NEEDED  IF  e  >  5XI0-13 


INCREASING 
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DERIVATIONS  OF: 

(18),  (19),  (20);  (48),  (57);  (59),  (60),  (61);  (69),  (70) 

In  this  appendix  derivations  are  given  for  equations  (18), 
(19),  (20);  (48),  (57);  (59),  (60),  (61);  (69),  (70). 

A.  Derivation  of  (18),  (19),  (20) 

Equation  (18)  is  given  by 

Ix(a,b)  =  x  Ix(a  -  l,b)  +  (1  -  x)  Ix(a,b  -  l)  .  (18) 

This  equation  is  proved  by  first  establishing  the  relation 


Ix(a  +  l,b  -  1)  =  Ix(a,b)  -  .  lia  +  b) 

r(a  +  i)  r(b) 

An  integration  by  parts  on  Bx(a,b)  gives 


xa  (l  -  x)b-1  . 


(107) 


Bx(a,b)  =  Bx(a  -  l,b  +  l)  -  f  "1 A  T .  (108) 


Therefore 


Bx(  a  +  l,b  -  1)  =  Bx(  a,b)  -  —  A  "/)  —  .  (109) 


b  -  1 


b  -  1 


Multiplying  (109)  by  r(a  +  b)/[r(a  +  l)  r(b  -  l)]  leads  directly 
to  (107). 

The  proof  for  (18)  follows.  From  (107) 


,a-  1 


Ix(a  -  l,b)  -  Ix(a,b  -  1)  =  rA  -  yr-1  r(a  +  b  -  l) 

r(a)  r(b) 


=  X-lA.A  -  x)b~.1  r(a  +  b) 
r(a)  r(b) 


1  - 


a  -  1 


b  -  1 


a  +  b-1  a  +  b-1 
Now  assuming  x  ^  0,1,  (110)  .may  be  written  as 


(110) 
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—  [ I  (a,b) ]  =  JL  [x  I  (a  -  l,b)] +  —  [(l  -  x)  I  (a,b  -  l) ]  . 
dx  dx  dx 

(111) 

Carrying  out  the  obvious  integration,  gives  (18)  plus  an  integra¬ 
tion  constant  which  can  be  shown  to  vanish  by  letting  x  tend  to 
zero.  It  is  obvious  that  (18)  also  holds  for  x  =  0  and  x  =  1, 
and  the  proof  is  complete. 

In  order  to  derive  (19),  the  following  relation  is  used: 
b  Ix(a,b  +  1)  +  a  Ix(a  +  l,b)  =  (a  +  b)  Ix(a,b).  (112) 
Equation  (112)  is  proved  by  writing  Bx(a,b)  as 

X 

Bx(a,b)  =  f  ta_1  (1  -  t)b_1  [(l  -  t)  +  t]dt 


=  Bx(  a,b  +  1)  +  Bx(a  +  l,b)  .  (113) 

Multiplying  (113)  by  r(a  +  b  +  l)/[r(a)  r(b)]  and  using  (45) 
gives  (112).  The  index  b  is  reduced  by  unity  throughout  (112); 
this  does  not  affect  the  validity  of  the  relation,  and  subsequently 
Ix(a,b)  as  given  by  (107)  is  substituted  for  the  second  term  on  the 
left  hand  side  of  (112).  The  result,,  after  some  trivial  algebra, 
is 


Ix(a,b)  =  Ix(a,b  -  1)  +  xa  d  "  x)b_1  >  *  1.  (114) 

When  b  =  1,  lx(a,0)  is  to  be  interpreted  as  zero. 

By  applying  (5)  to  (114),  or  by  manipulations  similar  to  those 
used  for  deriving  (114),  another  useful  result  is  obtained, 

I  (a,b)  =  I  (a  -  l,b)  -  Lif-tJB.  -1!  x3-1  (1  -  x)b  ,  a  >  1,  (115) 
x  x  r(a)  r(b) 

where  lx(0,b)  =0  in  (115).  Equations  (19),  (20)  follow  from  (107) 
by  setting  b  =  3/2  for  ( 19), and  by  setting  a  =  l/2  and  increasing 
b  by  unity  throughout  (107)  for  (20).  If  one  subtracts  (115)  from 
(114)  the  result  is  equivalent  to  (107). 
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B.  Derivation  of  (48),  (57) 

Equation  (48)  follows  easily  by  writing  (114)  as  a  tele¬ 
scoping  series,  where  b  is  replaced  by  a  running  index  i,  such  that 

3 

Ix(a,b)  =  Ix(a,j)  =  Z  [lx(a,i)  -  Ix(a,i  -  l)] 

i  =  i 

=  z  xa  £(„a  +  1  ~  ll  (1  -  x)i_1 
1=1  r(a)  r(i) 

which  is  also  (48). 

Equation  (57)  follows  easily  also  by  using  (45). 


ai+i  _  xa  (l  -  x)1  r(a  +  i)/[ r( a)  r(i  +  l)] 
al  xa  (l  -  x)1_1  r(a  +  i  -  l)/[r(a)  r(i)  ] 


=  (1  -  x)  .  (117) 


C.  Derivation  of  (59),  (60),  (6l) 

As  in  deriving  (48),  (59)  is  obtained  by  writing  (114) 
as  a  telescoping  sum.  However  in  this  case,  b  =  j  -  l/2  and 

j-1 

Ix(a,b)  -  Ix(a,l/2)  =  Z  [lx(a,i  +  l/2)  -  Ix(a,i  -  l/2)] 

1  =  1 


(116) 


Thus 


=  3Z  x-  ha-t  i  -  1/g?  (1 
1  =  1  r(  a)  r(i  +  1/2) 


(118) 


which  is  the  result  desired. 
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Equation  (601)  is  directly  deducible  from  (ll8)  by  setting 
a  =  l/2.  The  equation  (60a)  is  also  easily  proven  by  applying  (5) 
to  (60b).  The  term  Ix(l/2,  1/2)  is  given  by 

x 

i,(i/2,  1/2)  =  (r(i)/[r(i/2)  K1/2)]}  f  — M —  .  (119) 

J  /t  /n 

The  transformation 

t  =  sin2  0  y  0  s  0  <n/2  j 

applied  to  the  integral  of  (119)  gives  for  Ix(l/2,  l/2) 7 

sin"”1  Tx 

Ix(l/2,  1/2)  =  i  2  r  .£in_0  £2§JL  d@  =  2.  sin-1  ,/x  ,  (l20) 

J  sin  0  cos  0  ~ 

which  is  equivalent  to  (6l). 

D.  Derivation  of  (69),  (70) 

These  results  are  obtained  in  exactly  the  same  way  that 
(117)  was  generated. 
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DERIVATIONS  OF:  (91),  (92). 

PROOF  THAT  Un^r(u)  <:  Un^r(o)  ,  0  <;  u  <:  1 

In  this  appendix  the  three  following  results  are  proved: 

(A)  Equations  (9l),  (92) 

(b)  The  maxima  of  |un  r|  decrease  monotonically  as  a  function 
of  u  on  the  interval  [0,1]. 

A.  Proof  of  Equations  (91)  and  (92) 

In  (89),  Un  r  is  defined  accordingly 

un  r  =  SL  [(1  -  u2)n]  ,  0  £  u  «£  1  .  (89) 

’  du 

Equation  (94)  states  that  Un^r  is  given  hy 

[r/2] 

Un  r  =  Z  (-  l)r-1  — 2^-al_r.Ir(n  +  IV  ur“2i  (l  _  u2)n- 

3  i_0  il(r  -  2i) i  r(n  -  r  +  i  +  1) 


n  >  r  ,  (9l) 

where  [r/2]  represents  the  greatest  integer  in  r/2.  The  proof  is 
hy  induction.  Thus  for  r  =  0,  1,  2,(9l)  is  easily  seen  to  he 
valid.  It  is  necessary  to  show  (91)  holds  for  x  assuming  it 

holds  for  Un^r;  Un^r+1  would  he  given  hy 
fr+ll 

u  =  (_  !)r+l-i  2r~21+1  (r  +1)1  r(n  +  1) 

n,r+1  y  il(r  -  2i  +  l)l'  r(n  -  r  +  i) 

.  ur-2i+1  (1  -  u2)"-^1-1  .  (121) 

The  proof  follows: 

Introduce  ,  such  that 


=  (-  l)r-1  _ 2r~21  rlp(n  +  1) 

i!(r  -  2i)lp(n  -  r  +  i  +  l) 


i  ^  0;  Aj  =  0,  i  <  0 . 

(122) 


*r+l 

y 
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Now  differentiating  (9l),  and  subsequently  using  (122)  one  obtains 


A  [U  1  =  (-  l)r+1  g^rjntl)  ur+1  (1  -  u2)11-1-1 
du  n’ r  r(n  -  r) 

[r/a] 

+  L  [ ( r  -  2i  +  2)  A  -  2(n  -  r  +  i)  A  ]  ur_21+1 
i=l  1-1  1 

•  (l  -  ■u2)n~r+1”1  +  (r  -  2  [r/2])  Aj-^y^j  ur"2  fr/2^ 

•(1  -  u2)n-(r-[r/a])  .  (123) 

The  constant  factor  under  the  summation  sign  of  (123)  can  be 
simplified.  Thus 

(r  -  2i  +  2)  A^^  -  2(n  -  r  +  i)  At 

is  equal  to 

(-  i)r“1+1  _ gr-51+2  rlrfn  +  l) _ 

(i  -  l)l(r  -  2i  +  l)ip(n  -  r  +  i) 


+  (-  l)r“1+1  2r~gl+1  rlr(n  +  1) _ 

ii(r  -  2i)if(n  -  r  +  i) 

_  (_  p) r-i+i  2r~gl+1  (r  +  l)  I  r(  n  +  l) 
i!(r  -  2i  +  l)if(n  -  r  +  i) 


(124) 


The  last  term  in  (123)  is  equal  to  the 


odd.  The 


r  +  1 


r  +  1 


term  of  (l2l)  for  r 


term  of  (l2l)  for  r  even  is  included  in  the  [r/2] 


term  of  the  sum  in  (123).  It  therefore  follows  that  (123)  is  equal 
to  (121),  and  the  proof  is  complete. 
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=  -  2(n  +  1)  jr+1-  [u(l  -  u2)n] 
dur+1 

=  -  2(n  +  1)  [u  U^r  +  (r  +  1)  Un^r]  ,  (127) 

where  Leibnitz*  s  rule  was  employed  again  to  obtain  the  last  equa¬ 
tion.  Subtracting  both  sides  of  (127)  from  both  sides  of  (126) 
gives 

(1  -  u2)  U^r  +  2(n  -  r  -  1)  u  U^p  +  (r  +  l)  (2n  -  r)  Un#r  =  0, 

(92) 

which  is  equation  (92). 

B.  The  absolute  values  of  the  extrema  of  (89)  decrease  as 
a  function  of  u  on  [0,li  provided  n  >  r  +  l/2.  The  proof  for  this 
statement  was  suggested  by  techniques  used  by  Szego  in  [13,  Chap¬ 
ter  VII] . 

Consider  a  functon  f  such  that 

f  =  A  y^(u)  +  cp(u)(y'f  >  (128) 
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where 


A  is  a  positive  constant,  and  cp(u)  is  non-negative  for  u  in 
[0,1],  Therefore 


f  20, 


and 


f '  =  y'  [ 2cpy"  +  cpV  +  2 Ay]  .  (129) 


Now  let 

y  =  Un,r  ,  cp(u)  -  |  (1  -  u2)  ,  2A  =  (r  +  1)  (2n  -  r)  . 

(130) 


Then,  after  substituting  (92)  into  (129) 

f'=2y/[-(n-r-l/2)uy/],  0<u^l.  (l3l) 

Therefore  it  is  concluded 

f '  <  0  if  n>r  +  l/2,  u  ^  0,  y'  ^  0  ,  (132) 


where  the  inequality  of  (132)  also  insures  A  as  defined  in  (130) 
to  he  positive.  Since  f '  is  negative  on  (0,1)  at  points  where 
y'  i1  0,  this  means  that  f  is  a  decreasing  function  of  u  on  [0,1], 
The  clinching  argument  follows  by  considering  those  values  of  u 
for  which  y'(u)  =  0,  on  [0,1],  i.e.,  the  extrema  points  of  y, 
which  we  call  um.  For  such  u,  (128)  can  be  written  as 

^(Um)  =  f(um)/A  .  (133) 

But  since  f  is  a  decreasing  function  of  u,  then  y^u^  cannot 
increase  as  the  um  increase  from  0  to  1.  If  r  is  even,  u  =  0  is 
a  point  of  the  set  {um},  because  it  is  evident  from  (91)  that 
y'(0)  =  Un  r+  (0)  =  0.  If  r  is  odd,  u  =  0  does  not  belong  to  the 
set  {um}  .  Thus 

K,r(0)l  >  |un,  r(u)|,  0  <  u  s  1,  r  even  .  (134) 
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The  result  which  has  just  been  proved  can  also  tie  deduced  directly 
from  a  theorem  given  by  Tricomi,  [16;  p  99].  The  theorem  essentially 
states  that  if  a  differential  equation  has  the  form 


d 

du 


p(u)  $L 


+  P(u)  y  =  0  , 


(135) 


such  that 

a)  p(u)  and  P(u)  and  their  first  derivatives  are  continuous  on 
(a,  b),  i.e.,  p(u),  P(u)  are  in  C'  on  (a,  b), 

b)  [p(u)P(u)]  is  a  non-decreasing  (non- increasing)  function  of 
u  in  (a,b), 

c)  P(u)  /  0  in  (a,  b), 

then  the  maxima  and  minima  which  occur  in  (a,  b)  of  any  integral 
y(u)  of  (135)  are  such  that  the  corresponding  values  of  I y |  form  a 
non- increasing  (non-decreasing)  sequence.  If  the  hypotheses  of  this 
theorem  are  satisfied  on  the  half-open  interval  [a,  b),  then  it  is 
easily  shown  by  going  through  Tricomi's  proof  step  by  step  that  the 
conclusion  holds  on  the  half-open  interval,  that  is,  the  extrema  on 
[a,  b)  are  such  that  the  corresponding  values  of  |y]  form  a  non¬ 
increasing  (non-decreasing)  sequence. 


Equation  (92)  is  easily  put  in  the  form  of  (135),  (see  [16; 
p  96]),  so  that  (92)  becomes 


d_ 

du 


-(n-r-l) 

(1  -  u2) 

duJ 


-( n-r) 

+  (r  +  l)  (2n  -  r)  (l  -  u2)  y  =  0, 


where 

-  -(n-r-l ) 

p(u)  s  (l  -  u2)  ,  P(u)  =  (r  +  1)  (2n  -  r)  (l  -  u2) 


(136) 

-(n-r) 


(137) 

On  to,  1),  p(u)  and  P(u)  are  obviously  in  C P(u)  /  0,  and 
[p  (u)P(u)]  is  non-decreasing  provided  n  s  r  +1/2.  Therefore 
the  hypotheses  of  the  .modified  theorem  (for  the  half-open 
interval)  are  satisfied,  and  the  conclusion  of  the  modified 
theorem  holds  and  implies  the  result  which  was  to  be  proved. 
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